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ABSTRACT 


This investigation deals with an analytical study of the effects of 
hydrogen injection and chemical reaction on the flow properties of Couette 
flow with s|)ecial emphasis given the diffusion model assumed for the 
calculations. The two diffusion models chosen for the present analysis 
are Fick's law diffusion and multicomponent diffusion. For most boundary 
layer and Couette flow analyses the approximate Fick’s law diffusion model 
is used as it results in considerable mathematical and numerical simpli- 
fication over the more exact but cumbersome multicomponent diffusion 
model. There is discussed in the literature the use of Couette flow to 
simulate the two-dimensional laminar boundary layer, however, there is 
no literature available concerning hydrogen injection into chemically 
reacting Couette flow with property variations nor is there literature 
available on the effect of the diffusion models used for the Couette 
flow solutions. The purpose of this study was to obtain solutions for 
the chemically reacting Couette flow with variable transport properties 
and hydrogen injection for the two diffusion models. 

Solutions to the governing equations for^Couette flow were obtained 

■’4 
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^or the two diffusion models over a range of hydrogen injection rates. 

The results indicate that there are significant differences between the 



solutions for the two diffusion models and these differences are 
manifested most in the concentration profiles and the lower wall 


heating rates. 
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I. INTRODUCTION 

The reduction of intense aerodynamic heating encountered by reentry 
and hypersonic flight vehicles through the use of mass transfer cooling 
has become widely accepted. Whether this mass transfer cooling is 
accomplished by ablation or transpiration, the gases injected into the 
boundary layer are generally quite different from the main stream flow. 
Since the heating reduction is greatest with low molecular weight gases, 
molecular hydrogen is usually a major component of the injected gases 
especially in the ablation of polymeric materials. The introduction of 
hydrogen into boundary- layer flow complicates the analysis as large 
property variations occur and molecular diffusion and chemical reactions 
must be considered. 

In most analyses Pick's law diffusion is assumed as it is an easily 
applied approximation to the more exact but mathematically cumbersome 
multicomponent diffusion model. However, since the Pick's law diffusion 
model is an approximation the calculated diffusion velocities will be 
in error, especially when there are large differences in molecular weight 
of the diffusing species as is the case when hydrogen is present in an 
airstream. Thus a comparison of the diffusion models is necessary to 
provide an estimate of the errors incurred when using the approximate 
model. 

In making a comparison of the diffusion models any simplification 
that can be used without concealing the important aspects of hydrogen 
injection into an air boxmdary layer is desirable. In the literature 
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the one -dimensional Couette flow model has been used to simulate the 
two-dimensional laminar boundary layer. However, the sources available 
consider only hydrogen injection into an air Couette flow with constant 
properties and no chemical reactions. 

The present study is twofold in purpose: first the effects of 

hydrogen injection with property variations and chemical reactions will 
be considered and secondly comparisons will be made between two diffusion 
models. As in the literature, the one-dimensional Couette flow model will 
be used to simulate the two-dimensional laminar boundary layer. In this 
Couette flow representation the velocity of the moving plate represents 
the free- stream velocity, while the distance between the plates simulates 
the boundary -layer thickness. 

This investigation was conducted under the auspices of the 
National Aeronautics and Space Administration at the Langley Research 
Center, Hampton, Virginia. 


II. REVIEW OF LITERATURE 


There exists surprisingly little information in the literature 
concerning the effects of the diffusion model on the solutions obtained 
for a chemically reacting airflow with hydrogen injection. There are 
no direct comparisons between the approximate Pick's law diffusion model 
aiid the more exact multicomponent diffusion model available from the 
literature, however, the analysis of Libby and Pierucci (ref. l) does 
consider hydrogen injection into a laminar boundary layer with variable 
properties, chemical reactions and multicomponent diffusion, but these 
solutions are compared to rather limited constant property (Prandtl and 
Schmidt numbers equal to one) solutions, making the comparisons somewhat 
unrealistic and giving no insight into the effect of the diffusion model 
utilized. This thesis differs from the above analysis in several respects. 
First, the solutions for the approximate diffusion model analysis will 
employ the same assumptions as the multicomponent diffusion analysis, 
except for the diffusion model itself; and secondly, the present analysis 
uses the Couette flow model to simulate the two-dimensional laminar 
boundary layer. 

There are several sources of information in the literature on 
hydrogen injection into Couette flow with the principal analysis being 
that of Eckert and Schneider (ref. 2), but because of the assumptions 
of no chemical reactions and constant properties their solutions are of 
only limited usefulness. A variable property analysis is given by 
Simon, et al. in reference 3 where liydrogen is injected into an inert 
stream with no chemical reactions. The present analysis differs from 
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these latter two references in that variable properties, chemical 
reactions, and two diffusion models are considered. Also, the present 
analysis will not employ the flame sheet approximation as did Libby and 
Pierucci to define combustion but instead a diffusion flame will result 


from the solution of the governing equations. 






III. ANALYSIS 

Figure 1 shows the Couette flow model. The lower porous surface, 
at y = 0, is stationary while the upper porous surface, at y = s, moves 
with a uniform velocity u„. The lower surface is at the temperature 
Tw and the upper surface at Too* The hydrogen gas, initially at temper- 
ature Tw> is injected perpendicularly into the flow uniformly through 
the stationary surface, and removed uniformly through the upper surface. 

Basic Equations of Motion 

The basic equations governing motion in a multicomponent chemically 
reacting gas are taken from Scala, reference 4. These equations are the 
continuity equation. 


+ pV • V = 0 
Dt 


and the momentum equation^ 




Pi^bi 


i 


where it is the pressure tensor, 
and the energy equation. 


Dh 
P — 
Dt 


1 


( 1 ) 


( 2 ) 


( 5 ) 
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In this equation Q represents the energy flux due to temperature and 
concentration gradients anc. is given by 



The viscous dissipation function in equation (3) is given by the 

following relationship: 



and finally there is the species continuity equation 

DKi 

p 

Assumptions for Present Analysis 
For the present analysis the following assumptions are made: 

1. The flow is steady. 

2. The model is one dimensional. 

3 . Thermal diffusion effects are neglected. 

Diffusion stress effects are neglected. 

5* No body forces considered. 
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6. No pressure gradients. 

7- No radiation heat transfer. 

8. Chemical equilibrium exists throughout the flow. 

9* Gas properties depend on local mixture concentration and 
temperature . 

10. Both surfaces are impermeable to the main stream gas. 

11. Prandtl's order of magnitude analysis is applicable (ref. 

Equations of Motion for Couette Flow 
Using the above assumptions the basic governing equations of 
can be reduced to the following forms for Couette flow. 

Continuity equation, 

= 0 

ay 


momentum equation, 


du _ d 
dy dy 



energy equation, 


dh 

pv — = 
dy 



species continuity equation. 


dKi A 


5). 

motion 


(5) 


( 6 ) 


(7) 


( 8 ) 
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These equations are similar to the two-dimensional laminar boundary 
la yer_ equatlons except for the Couette flow being one dimensional. 

A further simplification of the energy equation can be obtained 
through the Introduction of: 


=I Vi 




Cpj^ dT + 


( 9 ) 


Thus 


dy 


I 


KfCpi 



dKj 

dy 


( 10 ) 


or 


dy 



(11) 


where 



The energy equation becomes: 


pvC 


^ A. 

PR dy ” dy 




PlVihi 


( 12 ) 
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A simplification of the species continuity equations can be obtained 
through the introduction of the concept of "elemental" mass fractions 
as expressed by Lees in reference 6. The "elemental" mass fraction 
concept results from the fact that the mass of individual chemical 
elements is preserved in any chemical reaction not involving nuclear 
transformation. The "elemental" mass fraction is given by the 
expression: 



The species continuity equations for the elements can be obtained by 
multiplying equations (8) by M 

results the "elemental" species equations 

The introduction of the "elemental" mass fraction eliminates the species 
production terms (cuj_) of equations (8) and reduces the number of 
calculations to be made. There is now one equation of the above form 
for each element as opposed to one equation of the form of equation (8) 
for each chemical species. 

In the present analysis there will be three elements H, N, and 0, 
and four chemical species ©2^ ^ 2 ^ 1 ^ 2 * ^2® considered with one chemical 

reaction of the form: 


summing over i, and there 






I- ^ 
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Hg + I Og HgO 

This same chemical system was used by Libby and Pierucci in reference 1, 
eind while it does not consider dissociation or ionization, the sx>ecies 
considered do represent the major products of hydrogen combustion in an 
airstream. Also, the species considered have the necessary variation in 
molecular weight which is essential to the diffusion model comparisons. 

Boiindary Conditions 

At the moving surface, (y = s) , the following boundary conditions 
apply: 

T = t;, 

Ki = Ki 

00 

In order to simulate the two-dimensional boundary layer the 
"elemental" mass fraction for hydrogen must be very small. This creates 
a correspondingly small "elemental" hydrogen density and since the 
continuity equation must be satisfied (pv = constant) the transverse 
velocity becomes very large. This introduces some uncertainty as the 
transverse velocity was assumed small in comparison with the main flow 
velocity in order to accomplish the reduction of the general equations 
of motion. This apparent contradiction is inherent in the use of the one 
dimensional Couette flow to simulate the two-dimensional laminar boundary 
layer and arises primarily from the assumption of the porous surfaces 
being permeable only to the hydrogen. 


LI 

At the lower surface, (y = O), the boundary conditions are: 

u = 0 
T 

The boundary conditions on the '‘elemental" species continuity equations are 
derived as follows. Integration of the global continuity equation yields 

pv = constant = (pv)^ (l4) 

Using this relation the "elemental" species continuity equation 
can be integrated to give: 


pv Kj 



Mj 

% 


PiVi = constant 


(15) 


The following subscript notation is adopted for the "elemental" species. 

element subscript 


0 

H 

N 


1 

2 

k 


Assuming Pick's law diffusion for illustration the diffusion velocities are 


PiVi = - pD 


dKi 

dy 


(15a) 


and assuming the same relation holds for the "elemental" species 


^ ^ dK^ 

p. V4 = - pD 

dy 


(15b) 
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Using the Fick's law relation, equation (l5a), the summation term of 
equation ( 15 ) becomes, 


V . \' M. dKi 

1 ki kii^^ 


(15c) 


The definition of the "elemental" mass fraction is 


~ V 

• ^ 


(I5d) 


Differentiating equation (ijd) one obtains. 


dK 


dy 


=y ^ 

Z_J M-I 


dKi 
Mi dy 


(I5e) 


Substituting equation (l5c) into equation (l5c) the following relation- 
ship is obtained. 


z 


dKi 

5lJ PiVi = - PD ^ 


or from equation (l5b). 


z 


M., 

^ij sf Pi^i 




Thus equation (15) becomes 




pvKj + pjVj = constant 


(I5f) 
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Considering first the injected species, hydrogen, equation (I5f) 
becomes; 

pgV + P2^2 “ constant 
P 2 (v + V 2 ) “ constant 

pgVg = constant (16) 

Evaluating the constant at the lower wall (y = O) the above equation 
becomes: 

P2V2 = (P2V2)^ = (pv)„ 

Thus the boundary condition on the "elemental" hydrogen continuity 
equation becomes; 


^ PlVl„ = (pv)„ 

or 




Kaw =.i- 


(pv)„VL Mi 


P^V 


i’i 


( 17 ) 


w 


A similar procedure is followed in evaluating the constant for the 
main streeim components where 


PlVi = (PiVi)^ 


= 0 


P4V4 = {puykh = 0 



Ik 


The boundary conditions for these elements 


K 


Iw 




(18) 


K 


4w = ■ 


' i ' w 


Nondimensional Form of the Governing Equations 
The following new variables are introduced: 


( 19 ) 


T) = 


1 

S 


U = 


e = 


*•00 


^ . pvs _ 

Oq 

M«) Moo 


The governing equations in nondimensional form are: 


momentum equation, 


6 ^ iL. ^ 

° <in 'in U» 'iny 


( 20 ) 


energy equation, 


'PR de d 


A d0\ ^ ^ ^ 


dU^ 


° Cpj^ dq Cpj^TeoM« \dTl 


-dl 


PiVis 




( 21 ) 
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species continuity equation, 


Bo 





piXif\ 

Mi / 


(22) 


Nondimensional Bounciary Conditions 


At 11=1 


U = 1 

e = 1 


Ki = Ki^ 


at q = 0 


U = 0 
0=9 

w 



Wy M,i PiYis\ 

Mi M.. / 


Klw = 



ii7 


Mj PiViS> 


Moo 


w 





t 



(25) 


(24) 


(25) 


Heat Treuisfer at Lower Wall 



(26) 
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Transforming equation (26), 


% = 


Qws 






PiVjs 

Mto 




w 


( 26 a) 


Shear Stress at Lower Wall 


du\ 

^ d^) 


Transforming aquation (27), 





MooMqo 


Moo dT] y 


(27) 


(27a) 


Solid Wall Couette Flow 

For the solid wall case, the mass transfer is zero. The continuity- 
equation integrates to: 

pv = 0 ( 28 ) 

The governing equations then reduce to; 
momentum equation, 

dy \ dy 
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energy equation, 


dy 



(30) 


On introduction of the nondimen sional variables these equations become; 
momentum equation, 



( 31 ) 


( 32 ) 


Integrating each equauion once and evaluating the constant at the 
lower surface the equi'.^ cions in nondimensional form are: 
momentum equation. 


dll [1 


( 33 ) 


where 




w 


energy equation, 

— \ - 1 M(T))Uoo/ dU\ ^ ^ ^2 

dq j Mq) Jq \dq / Mq) 


(3^) 
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where 



Boundary Conditions 


T) = 1 


U = 1 

e = 1 


n = 0 


u = 0 
e = 


Heat Transfer at Lower Wall 


dT 


w 


(55) 


Transforming, 


QwS _ \i d9\ 


(35a) 


Shear Stress at Lower Wall 


du\ 


(56) 
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Transforming, 



dU\ 

MooMoo Moo dr\) 

w 


(36a) 



IV. GAS PROPERTIES 


The thermodynamic and transport properties are calculated by the 
met nods listed below. The gas mixture is assumed to be at one atmosphere 
prctsure for all calculations. 


Chemical Composition 

The following reaction is considered for the present analysis: 


H2 



HgO 


In addition, N2 is present in the main Couette flow giving a total of 
four chemical species to be considered in the equilibrium calculations. 
The chemical equilibrium equations will be formulated in terms of the 
"elemental" mass fractions to facilitate calculation of the equilibrium 
composition from the solutions of the species continuity equations. 



(37) 

(38) 


(39) 
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Dividing equation (^) into equation ( 58 ) and rearranging one obtains the 
following equation for the mole fraction of N 2 : 



^H20 

2A 




where: 


% 

Mo% 


Dividing equation (*(0) into equation (39) anii rearranging one obtains: 




(42) 


Where: 


Mjj % 

Substituting equation (i*-l) into equation (^+2) the equation for the mole 
fraction of H 2 is obtained: 



(^ 3 ) 


C 


B 

A 


where : 
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Using the relation 


^ X, = 1, 


Substituting equations (43) and (4l) into equation (44) the following 
relationship results: 


( 44 ) 




V = 1 


(**5) 


Letting 


-{■ 


E = U + C + 


= 

\^2 2AJ 


■G =i 


X„ 0 = (G - E®0,) 


(‘* 6 ) 


Combining equation (46) and (43) one obtains: 


° " “°2 


(‘^7) 


The equilibriiom constant is related to the mole fractions by 


pl/2 


S = 


%o0 


XH2(XQ2)^'^^ 


(48) 
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Substituting equations (^) and (47) into equation (48) there results 
the following relation: 

(G - EGXog) 

{(£ - l)(G - EGXo^) + CX02] (Xq2^/2) 



Letting 


I = P^/% 

P 



the above relation reduces to 


I = 


(G - BGXog) 


• (HG - HEGXog + CXo2)(Xo2^/2) 


IHGXo,,^/^ - = G - EGXq, 


(- - eg)xq 5/2 + 52 xq + gxqJ-/^ = — 

\H /^2 IH ^2 ^2 IH 


(*^9) 


Since all the constants (A through I) in equation (49) are known 
quantities, equation (49) can be solved iteratively for the unknown mole 
fraction Xq^* With the mole fraction Xq^ the other mole fractions can 
be easily determined from equations (46), (43), and (4l) . The eqailibrium 



constant used in these calculations is taken from the JANAF tables, 
reference 1, where it is tabulated in one hundred degree Kelvin increments. 
The molecular weight of the mixture is determined from: 


M =y 


( 50 ) 


The mass fractions are determined from the mole fractions by: 



XjMi 

M 


( 51 ) 


Thermodynamic Properties 

The mixture density is obtained from the perfect gas law 


P 


m 

RT 


( 52 ) 


and the enthalpy of the individual species is taken from the JANAF tables, 
reference 'J , and the mixture enthalpy is calculated by 



( 53 ) 


where 




dT + h^^ (ref. 7) 


The reference temperature for the present calculations is 0° K. 
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The specific heat of the individual species is also taken from 
reference 7 . The specific heat of the mixture is obtained from the 
relation 



( 5 ^) 


The derivative of Ki is found numerically by solving for the K^'s at 
a temperature 25^ K above and below the given temperature with: 


dKi 

dT 50 


(55) 


A comparison of the present method and the method of Zeleznik and 
Grordon, reference 8, is given in table I which shows very good agreement 
between the two methods. 


Transport Propierties 

Rigorous kinetic theory expressions for the viscosity and thermal 

conductivity of gas mixtures have been developed and are presented by 
Hirschfelder et al. in reference 9^ t)ut these expressions are mathemati- 
cally cumbersome. Somewhat simpler relations, which are approximations 
derived from the rigouous expressions are given by Brokaw in reference 10 
and are used in the present analysis. These approximations are very accu- 
rate at low temperatures but the accuracy is expected to deteriorate some- 
what at higher temperatures due to uncertainties in the approximations 
and molecular constants used, however, these approximations are of suffi- 
cient accuracy for use in the exi)ected temperature range of the present 


analysis. 
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Mixture Viscosity 

The mixture viscosity is calculated from the pure component 
viscosities with the relation 



( 56 ) 


The coefficients are a function of the pure component viscosities 

and molecular weight ratios 



\ \“i/ 


2\fi 1 + 


Mi\\l/2 


( 57 ) 


For use in the present analysis the pure component viscosities are taken 
from Svehla, reference 11. 

Mixture Thermal Conductivity 

The mixture thermal conductivity is obtained from the relation: 

(58) 

where is the treuisfer of energy due to the translational motion of 

the molecules and is the transfer of energy between internal degrees 

of freedom and translational motion in polyatomic molecules. The 


monatonic mixture conductivity is obtained from the pure component 
monatonic conductivities with the relation 



( 59 ) 


The coefficient is obtained from the viscosity coefficient 

by following relationship 


'I'lJ = ^ + 2.41 


- 0.142 M.) 


(Mi + Mj)' 


(60) 


The internal mixture conductivity is obtained from the pure component 
internal conductivities with the relation 



(61) 


The pure component thermal conductivities are obtained from reference 11. 

The pure component transport properties used in the analysis are 
calculated in reference 11 by: 

Viscosity 


^ ^ 26.693N/^^10.6 


a%i(2,2)* 


(62) 
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Monatonlc Thermal Conductivity 


X' = M (65) 

m 


Internal Thermal Conductivity 



In calculating the viscosities and thermal conductivities using the 
above equations the Lennard- Jones (l2-6) potential was used to obtain 
the reduced collision integrals fi(2,2)*. The reduced collision 
integrals are given in table form in reference 9* The molecular force 
constants (see table II) used in the Lennard- Jones (l2-6) potential were 
obtained from experimental viscosity and thermal conductivities where 
possible and estimated by Svehla in reference 11 using empirical 
relations where no data was available. In the calculation of the reduced 
collision integral the molecules were all assumed to be nonpolar. This 
assumption does not introduce appreciable error since the force constants 
were obtained from viscosity measurements sind since the dipole-dipole 
interaction is negligible in higher temperature high energy collisions. 

A more complete discussion can be found in reference 11. 

Diffusion Transport 

The purpose of the present analysis is to compare solutions to the 
governing equations for Couette flow using two different diffusion models; 
the approximate Pick's law diffusion model and the more exact multicomponent 
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diffusion model. As will be seen below, the more exact multicomponent 
diffusion model entails a considerable number of mathematical operations 


and from a numerical analysis standpoint is not as desirable as the 
simpler but approximate Fick's law model. These two diffusion models are 
presented in the following sections. 

Multicomponent Diffusion 

The multicomponent diffusion fluxes were calculated using the 
following relations from reference 9* 


v-1 



Equation (65) is the Stephan -Maxwell relationship for the multicomponent 
diffusion velocities. 

V 

y = 0 (66) 

i=l 


Equation (65) can be rearranged to a more convenient form: 



(67) 


Multiplying equation (67) by and introducing the nondimen sional 

coordinates we have: 


ax, 

=7 

Meo dT) ^ 


x.x. 

i J 




Moo 


PlViS X,Xj 


Moo 


I 


J=1 ^ 


(68) 
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Similarly multiplying equation (66) by 


I 

i=l 


PiV,^ 


= 0 


(69) 


For the v component gas mixture the diffusion fluxes, 

are obtained from the simultaneous solution of v - 1 relations of the 

form of equation ( 68 ) and the relation given by equation (69) • 

The binary diffusion coefficients are calculated using the following 
relation from reference 9* 


'ij 


= 0.002628 


1/2 


t5(Mj^ + Mj) 
p(oij)2 nij(i’i)* 


(70) 


Again the reduced collision Integral, ^^ij is based on the 

Lennard -Jones potential and is taken from reference 9* The molecular 
constants used for these calculations and to obtain the reduced collision 
integrals are given in table II. The binary diffusion coefficients 
obtained from the above relation are shown in figure 2. 

Fick's Law Diffusion 

The Fick’s law diffusion fluxes are calculated according to the 
following relation: 


dKi 

p.V. = - pD — i 
dy 


(71) 
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Knuth in reference 12 states that a sufficient condition for the 
applicability of equation ( 71 ) is that the binary diffusion coefficients 
are equal to each other and to the Pick’s law diffusion coefficient. 

This assumption makes the Pick's law diffusion coefficient a pseudo 
binary diffusion coefficient and in the literature Pick's law diffusion 
is generally referred to as binary diffusion because of the appearance 
of equation (71)* The term binary diffusion will be adopted here for 
discussion purposes. 

Multiplying equation (71) by l/u„ and introducing the nondimen sional 
coordinates equation (71) becomes 

PiViS pD dK^ 


The diffusion coefficient is considered to be the same for all species 
and as such is a self -diffusion coefficient given by the following 
relation from reference 9i 


D 


0.002628 


(t^/m)^/^ 


(72) 


where the molecular constants are an average of those in table II. 
Again the Lennard -Jones collision integral is used. The mixture 
molecular weight is used in place of an average molecular weight in 
order to more accurately represent the diffusion process. Although 
not plotted on figure 2, the average diffusion coefficient lies above 


the lower set of curves and below the upper set. Thus at a given 
temperature the average diffusion coefficient lies in a narrow band 
of values depending on the mixture molecular weight. 


V. COMPUTATION 


The philosophy on the numerical analysis was to keep calculations 
straightforward and as simple as possible. Some changes were made in 
the numerical technique during the debugging process but overall the 
numerical analysis is straightforward while not always simple. 

Solid Wall Couette Flow 

The governing equations for the solid wall Couette flow are arranged 
as follows; 
momentum equation, 

^ = fu(’l) (75) 

dT) 


where 


fu(n) 



energy equation, 



(7»*) 


where 




dT) + 


C2 
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The solution of each individual equation is obtained using the 
"corrector” method given by Hamming in reference 13. This method is 
based on an iterative finite-difference procedure using the following 
equation: 



where Aq is the distance between finite-difference atatlons. According 
to references 13 and 1^ the iterated corrector process always converges. 

In the evaluation of the f^Cn) and fgCr]) it is necessary to obtain 
derivatives of several functions and this numerical differentiation is 
accomplished through the use of the equations given in appendix A. 

Initial starting profiles are obtained from constant property solutions 
to the governing equations for Couette flow. Then solutions to the 
energy and momentum equations are obtained by repeated (iterated) 
application of the corrector equation until the following error criterion 
is met: 

I - f/| < 0.000001 

That is, the right-hand sides of equations (73) and (7^) are evaluated 
at each finite -difference station and the resulting functions integrated 
using equation (75) • This process is repeated until the above error 
criterion is met. The simultaneous solution of the energy and momentum 
equations is accomplished by the iteration technique of Smith and Clutter, 


35 




reference 15* This iteration process is shown in figure 3 and the block 
labeled "Fluid Properties" is riiaply the determination of the transport 
properties. The simultaneous solution is assumed to converge when the 
momentum equation has satisfied the following error criterion: 

1%^^^ - Un^I < 0.000001 

The present method for the solid wall Couette flow has been compared 
to several calculations from the literature. Air Couette flow 
comparisons have been made with the constant property analytical method 
of Eckert and Schneider, reference 2, for the following conditions: 

T = 2l8° K 
00 

Mach No. =12 

e„ = 6 

Only the temperature profile is compared as the velocity profile is 
linear in both methods. Good agreement was obtained between the two 
methods as shown in figure 4. In the constant property solutions the 
terms "edge" and "wall" refer to the moving and stationary surfaces, 
respectively, where the values of the properties were fixed for the 
solutions. Also present in figure 4 for general interest is the variable 
property solution and it is seen that the constant property solutions 
do not represent the variable property solution to any great degree. 

A variable property solution for nitrogen Couette flow was taken 
from reference 3 and is compared to the present method for a nitrogen 
stream in figure 5^ with good agreement again being obtained. The 



velocity profile (not shown) also showed good agreement. This comparison 
was made for the following conditions: 

= 218° K 
Mach No. =12 
Tw = 872 ° K 

The nondimen sional coordinate 7 is given by 

a 

’’ ’ % 


where 



and the nondimen sional temperature 9-j_ is given by 



Couette Flow With Injection 

The governing equations for Couette flow with hydrogen injection 
are arranged as follows: 


momentum equation, 
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where 


fu(n) = 


-L 5H' 

Bq dT^lnoo dTJ 


(77) 


energy equation, 


"i - ^e(n) 


(78) 


where 


f§(Tl) = 


A de 






^ y (79) 




species continuity equation 


^ = . aV ^ ^ 


dT) 


dT] Z_j 
i 


Mco 


(80) 


The "elemental" species continuity equation can be integrated to provide 
a more easily applied form so that the "elemental" mass fractions can 
be directly determined. 
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Integrating: 


r\j 


6o Kj = - 


I 


+ constant 

Mi Mtx, 


Evaluating the constant for the main stream components at t] = 0 


for nitrogen 


constant = = 0 


for oxygen 

constant = = 0 

The constant of zero for both "elements" results directly from the 
boundary conditions, equations (24) and (25). 

Thus tae species continuity equation for the main stream components 
becomes: 


~ . . ±(y ^ 

J 5 o\^ Mi M» / 


(81) 


Since there are three "elements" in the system only two species 
continuity equations need be solved, for the sum of the "elemental" 
mass fractions equal one. 

Because some difficulty was encountered in the solution of the 
governing equations, the presentation of the mamerical analysis will 
be somewhat in a reverse order from that of the solid wall case. However, 
the Initial starting profiles are again obtained from the constant 
property solutions to the governing equations for Couette flow. 
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The simultaneous solution of the governing equations is again 
accomplished by the iteration technique of Smith and Clutter, 
reference 15, but with one minor change. The iteration process is shown 
in figure 5 and the block labeled "Fluid Properties" is expanded in 
figure 6. The "Fluid Proi)erties" loop is necessary because the iterative 
solution of the "elemental" species continuity equations requires a 
new calculation of the chemical composition and diffusion velocities 
after each iteration. The "elemental" si>ecies continuity equations are 
assumed to be satisfied when the boundary conditions meet the following 
criterion: 


^CAL 




< 0.002 


The choice of 0.2 percent as the error criterion was dictated by the 
economics of computer usage for it was found that the major iteration 
(time-consuming calculations) were those in the "Fluid Properties" loop. 
The overall simultaneous solution of the governing equations was said 
to have been obtained when the following error criterion on the momentum 
equation was met; 


Ui 


I+l 


N 





< 0.001 


The error criterion of 0.1 percent was again dictated by the computer 
time used to obtain a solution. Within the major loops of the iteration 
process in figure 2 the energy equation was assumed converged when; 


A . : 




- - 


ho 


A I+l A I 

9/-1 


< 0.001 


and again the error criterion was 0.1 percent. 

Initially the solution of the energy, momentum, and species con- 
tinuity equations was attempted using the "corrector" equation given 
previously and the three-point derivative formulas. However, when the 
simultaneous solution converged there were present large oscillations 
in the solution, especially in the temperature (9) profile. An attempt 
to alleviate the problem was made by attempting to use a more accurate 
derivative formula, such as the five-point formula given below from 
reference l6. 


^ ' ^N+2 ^^N+1 ~ ^^N-1 ^N-2 .gg. 

dqj^ " 12AT1 

where the error term was of order (Aq)^. The oscillations became worse. 
Although the above formula is more accurate than the three-point formulas 
of appendix A both formulas have the same feature; they do not include 
the value at the point of calculation. In place of the three- and five- 
point formulas a four-point formula was derived, see appendix B, which 
considerably reduced the oscillations but did not completely eliminate 
them. The "corrector" formula, equation (75) ^ is but a two-point 
(trapezoidal) integration relation and in the region of rapid changes 
in the function to be integrated it is not very accurate. However, 
reference 17 gives some six-point integration formul is that are more 
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accurate and these formulas, see appendix C, are used for the present 
calculations resulting in a smoothing out of the oscillations in the 
converged solution. 

The solution to the individual momentum and energy equations was 
obtained by iterating the six-point integration formulas until the 
following error criterion was met: 

0.000001 

That is, the right-hand sides of equations (76) and (78) are evaluated 
at each finite-difference station and the resulting functions then 
integrated using the six-point formulas of appendix C to obtain the 
velocity (U) and temperature ( 0 ) profiles. This procedure is repeated 
until the velocity and temperature profiles meet the above error 
criterion. 

The present method for Couette flow with hydrogen injection has 
been compared with the solutions of Simon et al., reference for the 
following conditions: 


T„ = 2 l 8 ° K 
Mach No. =12 
T^ = 872° K 
Re^ =0.5 

where : 



The stream is all nitrogen and has variable properties and binary- 
diffusion. In general very good agreement between the methods was 
obtained. Figure 7 shows the temperature profile comparisons with 
the agreement not as good as expected but within reason. Better 
agreement was obtained between the velocity and concentration profile 
as shown in figure 8. On the basis of the comparison made here it 
is assumed that the numerical technique is sufficiently accurate to 
carry out the present investigation. 


VI. DISCUSSION OF RESULTS 


The following values of the independent variables were used to 
obtain both binary (Pick's law) and multicomponent diffusion solutions: 


T„ = 218° K 
Mach No. = 6 / 
T^ = 872 ° K 

J 


For all cases 


and 5 q =0., 0.05, 0.1, 0.15, 0.2, 0.35, 0 . 5 , 0.75, 1-, 1.5* 

In the numerical calculations fifty finite -difference stations were 
used for all cases except the no injection case where forty stations 
were used. The solutions were obtained on a CDC 660O computer and run 
times for a single case varied from a few seconds for the no injection 
solution to about thirty minutes for the multicomponent diffusion 
solutions at large injection rates. The convergence of the total 
solution was fairly rapid, requiring generally about four iterations of 
the momentum equation. Referring to figure 3, each iteration of the 
momentum equation required fewer iterations of the energy equation, with 
the total number of energy equation iterations being about five times the 
total number of momentum equation iterations. By far the most iterations 
were made in the "Fluid Properties" loop of figure 6 with a total number 
of iterations varying from several hundred to several thousand depending 
on the initial starting profiles and the magnitude of the injection rate. 
The numerical technique of this thesis is not very sophisticated in 
terms of present day numerical analysis but there has not been a single 
case of nonconvergence encountered with this technique. 
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The no injection temperature and velocity profiles are given in 
figure 9 and it should be noted that the stream temperature Increases 
only slightly above the wall value, indicating that for the present 
conditions the wall temperature is less than but close to the adiabatic 
wall temperature. The velocity profile is not a linear profile due to 
the viscosity variation through the stream. 

Velocity Profiles 

The nondimensional velocity profiles are given in figures 1.0 
through l4. There does not appear to be any major difference between 
the solutions for the binary and multicomponent diffusion models, 
especially at the lower injection rates where the amount of hydrogen 
and water are substantially reduced in comparison to the oxygen and 
nitrogen. Also, the binary and multicomponent diffusion model velocity 
profiles in the region of the lower wall do not show any significant 
differences and in most cases the differences are hardly detectable. 
However, there is a difference in the nondimensional shear stress at 
the lower wall as shown in figure 15 . The shear stress for the multi- 
component diffusion model is higher than the corresponding binary 
diffusion model solutions for all injection rates. This shear stress 
difference results primarily from the mixture viscosity variations 
between the two diffusion models. The pure component viscosities for 
hydrogen and water are lower than those for nitrogen and oxygen, result- 
ing in decreasing mixture viscosity with increasing hydrogen and water 
concentrations. As will be shown in the section on concentration 


profiles, the hydrogen and water concentrations at the lower wall for 


the binary diffusion model are greater than the corresponding multi- 
component diffusion model solutions, hence, the mixture viscosity is 
less at the wall in the binary diffusion solutions. This decreased 
viscosity causes the somewhat reduced shear stress for the binary 
diffusion model as seen in figure 15 * 

In the present solutions the shear stress for both diffusion models 
increases substantially over the no injection condition in the inter- 
mediate injection rate range. This effect was not seen in the results 
of Eckert and Schneider, reference 2, where the shear stress was shown 
to decrease from the no injection condition for all hydrogen injection 
rates. Since reference 2 did not consider variable properties nor 
chemical reactions the increased shear stress of the present calculations 
is attributed to the inclusion of these phenomena in the analysis. 

Temperature Profiles 

The nondimensional temperature profiles are given in figures l6 
through 20. The differences between the binary and multicomponent 
diffusion models are greater for the temperature profiles than was the 
case for the velocity profiles with these differences being greatest at 
the intermediate hydrogen injection rates. The overwhelming effect of 
chemical reaction is seen by comparing the no injection temperature 
profile of figure 9 with figures l6 through 20. The increase in peak 
stream temperature over the no injection solution approaches a factor 
of three at the higher injection rates. Referring to figure 21 it is 
apparent that at low and intermediate injection rates the lower wall 
heating substantially increases but the heating rate then decreases as 
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the injection rate is increased further. These heating rate curves 
point out one of the largest differences between the diffusion model 
assumptions. The heating rates for the binary diffusion model are 
always larger than the corresponding multicomponent diffusion model and 
are positive, whereas the multicomponent diffusion model lower wall 
heating rate is negative at the largest injection rate. 

The no chemical reaction results of Eckert and Schneider, 
reference 2, do not show the large temperature and heating rate Increases 
that result when the exothermic hydrogen oxygen reaction is considered. 
Aside from this large effect of chemical reaction on the temperature 
profiles, the differences between the two diffusion models is more 
apparent especially in the heat transfer rates where the binary diffusion 
model generally predicts much larger heating rates than the multicomponent 
diffusion model. 


Concentration Profiles 

The differences between the two diffusion models is best seen in the 
concentration profiles of figures 22 through 30, where not only are there 
differences in the relative amounts of species but there are also some 
profile shax>e variations. The biggest difference occurs in the hydrogen 
concentration profile with the wall concentration reflecting this 
difference the most. Figure 31 gives the hydrogen concentration at the 
wall for both diffusion models and it is readily seen that the binary 
diffusion model concentration is much larger than the corresponding 
multicomponent model. This large difference alters the mixture transport 
properties at the wall, for hydrogen has a larger thermal conductivity 



and lower viscosity than the other species. This alteration of transport 
properties is the primary cause of the heating rate differences of 
figure 21, for the temperature profiles of the .two diffusion models are 
very close at the wall whereas the heating rate due to conduction alone 
is less for the multicomponent model due to the lower hydrogen 
concentration . 

The differences between the concentration profiles for the two 
diffusion models is primarily due to the increased diffusion velocities 
of the multicomponent model which means that a smaller chemical species 
gradient is needed to produce the same diffusion velocity as the binary 
model. This lowers the hydrogen concentration at the wall and also 
causes the reaction zone to be at a greater distance from the lower wall. 
This is best seen in figure 29 for an injection rate 6 q = 0.1 where 
for the binary diffusion model the reaction zone is away from the surface. 

The concentration profiles point out the area of greatest difference 
between the two diffusion models, with the concentration of hydrogen at 
the lower wall best representing these differences. 


VII. CONCLUSIONS 


The foregoing analysis has resulted in the following conclusions: 

1. Hydrogen injection into air Couette flow with chemical reaction 
alters the stream properties and causes changes in the major flow 
variables. 

2. There are significant differences between the two diffusion 
models, and the approximate Pick’s law diffusion model should be avoided 
except where gross approximations to the flow variables are needed. 

3 . The injection of hydrogen into air Couette flow with chemical 
reactions generally results in an increase in skin friction and heating 
rate . 

4. Major differences between the two diffusion models is manifested 
in the concentration profiles where not only are the concentrations 
different but there is some variation in profile shapes. 

5 . The multicomponent diffusion model solutions for the lower 
wall shear stress show an increase in the shear stress over the 
corresponding Pick’s law diffusion solutions. 

6. The multicomponent diffusion model solutions for the lower wall 
heating rates generally show a decrease in heating rate over the 
corresponding Pick’s law diffusion solutions. 


VIII. RECOMMENDATION 


It is recommended that a study similar to this one be carried out 
for the two-dimensional laminar boundary layer to establish the validity 
of using Couette flow with hydrogen injection and chemical reactions to 
simulate a chemically reacting two-dimensional laminar boundary layer 
with hydrogen injection. 
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APPENDIX A 


THREE-POINT laiRIVATIVE FORMULAS 

Figure 52 shows the finite -difference stations in the Couette flow 
and for the present analysis these stations are evenly spaced. The 
distance between stations is At). The finite-difference form of the 
first derivative is obtained by a Taylor series expansion about station N 
evaluated at stations N+1 and N-1, with 



_ %+l " ^N-1 (Aq)^ / d^ \ 

2An ■ 6 

The first derivative is thus correct to texms of order (Arj)^. 

The first derivative at the lower boundary (N = l) is also 
determined by the Taylor series expansion about station N but is 
evaluated at stations N+1 and N+2, with 


%+l = % + ^ 


dP\ ^ (An)^ / d% \ ^ (An)^/d^ \ 

2 \dTi2/jj 6 \dT]5/jj 




and 



The first derivative is again correct to terms of order 

In a similar manner solving for the first derivative at the upper 
boundary (N = NT) : 


dn ) 
'N=NT 


+ Fjj,2 - 4Fn.i 

Srj 


1 1 (^)}(^\ 

wn V H=iiT 


Likewise the first derivative at the upper boundary is correct to terms 
of order (Aq)^. 

The above formulas are used in the solid wall calculations wherever 
a first derivative is needed. 


APPENDIX B 


FOUR-POINT DERIVATIVE FORMULAS 


As in the case of the three-point derivative formulas the finite- 
difference points are evenly spaced as shown in figure 32. The finite- 
difference form of the first derivative at station N is obtained by a 
Taylor series expansion about station N and is evaluated at stations 
N-1, N+1, N+2. Thus: 



N+2 


= F 


N 




+ - - - - 


Solving for the first derivative 


dT] u 




The first derivative is correct to terms of order (^)5. The above 
equation applies in the InterveLL 2 £ N £ NT - 2. 

At station N = MT - 1 the first derivative is again obtained from 
the Taylor series expansion but is evaluated at stations N-2, N-1, 
and N+1, giving: 
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Again the first derivative is correct to terms of order (/^)5 and 
the above equation applies only at N « IfT - 1. 

At the lower boundary (N = 1) the first derivative is again obtained 
from a Taylor series expansion about station 1 but evaluated at stations 
N+1, N+2, and N+5* Thus: 


dT]2 \dr]V^ o 


+ - _ - > 


\ Vn \ ' 'n ' 'n 
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Solving for the first derivative — 

Again the first derivative is correct to terms of order In a 

similar manner solving for the first derivative at the upper boundary 
(N = NT); 

S)„ ■ ■ "'“-5} * 

Likewise the first derivative at the upper boundary is correct to terms 
of order (At])5. 

The above equations are used in the calculation of a first derivative 
wherever one is needed in the solution of the equations governing Couette 
flow with injection. 



APPENDIX C 


SIX-POINT INTEGRATION FORMULAS 


The six-point integration formulas given by Milne in reference 17 
and rearranged for the present analysis are: 


N = 2 


= Ft + —^1 495. fi + 1537. fp - 618. fz + 502. fL 

^ l 440 . p H- 

- 82. f5 + 9. fgj 


where the error term is: 


-862. (An)'^/ d'^F\ 

6o,48o. \dT) ^ / 


N = 3 


F3 =F2 


+ ^ — <- 4050. fT + 65430. fp + 75780. fz, - 7020. fk 

129600. 


- 1170. f^ + 630. f5 


with the error being: 


-77^1. 

6o48. 
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From N = 4 to NT - 2 


''n = %-l 144^. ^N-5 " ^N-J 

+ 8020 fjj - 950* ^N+1 "** ^10* 


the error term is: 


i+2356. \dii'/ 


N + NT - 1 


N = %-i ^ n;So:{‘ 

+ 10220 . %_i + 6370. % - 270. %+;^ 


the error term is: 



N = NT 



25920. 





- 14364. fjj_2 + 25686. + 8550. 


+ 8020. %_! 


- 2580. %_2 


, + 8676. %.5 
n} 



6o 


< ' 

the error term is: 


^ ‘> 5027 - 
2(il60. 
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Figure 5*- Main iteration process. 
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Figure 4.- Comparison for constant property Couette flow 

with no injection. 





65 


X 





Figure 5*- Comparison for variable property Couette flow 

with no injection. 
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This iterative 
procedure 
repeated until 
the boundary 
conditions of 
the species 
continuity 
equations ore 
satisfied . 


Figure 6.- Fluid properties Iteration process. 
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Figure 7«- Temperature profile comparison for nitrogen Couette 

flow with hydrogen injection. 
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(a) 8^, = 1.3 



U 


(b) 8 q = 1.0 

Figure 10.- Velocity profiles for 6 q = 1*5 and 1.0. 






(b) 6o = 0.5 


Figure 11.- Velocity profiles for 


6o 



= 0.75 and 0.5* 






(b) &0 = 0.2 

Figure 12.- Velocity profiles for 6 q =0.35 aud 0.2 





Figure 13.- Velocity profiles for 6 q 


= 0.13 and 0 . 1 . 
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Figure 19*- Temperature profiles for Oq = 0.13 and 0.1 
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Figure 21 


Heating rate results. 
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(a) Miilti component diffusion 



(b) Binary diffusion 

Figure 29.- Concentration profiles for 8^ = 0.1. 
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TABLE I.- COMPARISON OF EQUILIBRIUM CHEMISTRY RESULTS 


Temp. 




*H2 


XH2O 


Xr2 


X02 


1000 

Ref. 8 

Present 

method 

123 

123 

0.2800 

0.2799 

0 

0 

0.02967 

0.02961 

0.77885 

0.77886 

0.19148 

0.19153 


2000 


Present 


5000 


4000 


Present 

method 

, 

Ref. 8 

Present 

method 

740.9 

1093.7 

1094.5 

0.3382 

0.3561 

0.3566 

0.00277 

0.01653 

0.01659 

0.02679 

0.01289 

0.01278 

0.77778 

0.77241 

0.77240 

0.19265 

0.19817 

i 0.19823 


Initial conditions 


h = Cal/gm 

p = 1 atm 

CpR = Cal/gm-^K 

Kq = 0.2318 

T = °K 

Kh = 0.0021 






















































































